arXiv: 1506.01518vl [astro-ph.GA] 4Jun2015 


Mon. Not. R. Astron. Soc. 000 , ^|9] (2015) Printed 5 June 2015 (MN style file v2.2) 

Star Formation in Herschel *’s Monsters versus 
Semi-Analytic Models 

C. GruppionFf, F. Calura 1 , F. Pozzi 2 , I. Delvecchio 2 , S. Berta 3 , G. De Lucia 4 , F. 
Fontanot 4 , A. Franceschini 5 , L. Marchetti 6 , N. Menci 7 , P. Monaco 8 , M. Vaccari 9 

1 Istituto Nazionale di Astrofisica - Osservatorio Astronomico di Bologna, via Ranzani 1, 1-40127 Bologna, Italy. 

2 Dipartimento di Fisica e Astronomia, Universit.d di Bologna, viale Berti Pichat 6, 1-40127 Bologna, Italy. 

3 Max-Planck-Institut fur Extraterrestrische Physik (MPE), Postfach 1312, D-85741 Garching, Germany. 

4 Istituto Nazionale di Astrofisica - Osservatorio Astronomico di Trieste, Via Tiepolo 11, 1-34143 Trieste, Italy. 

5 Dipartimento di Astronomia, Universita di Padova, vicolo dell’Osservatorio 3, 1-35122 Padova, Italy. 

6 Department of Physical Sciences, The Open University, Milton Keynes MK7 6A A, UK. 

' INAF Osservatorio Astronomico di Roma, via di Frascati 33, 1-00040 Monte Porzio Catone, Italy. 

8 Dipartimento di Fisica, Sezione di Astronomia, via Tiepolo 11, 1-34143 Trieste, Italy. 

9 Astrophysics Group, Department of Physics, University of Western Cape, Bellville 7535, Cape Town, South Africa. 


Accepted 2015, May 27. Received 2015, May 6; in original form 2015, March 18. 


ABSTRACT 

We present a direct comparison between the observed star formation rate functions 
(SFRF) and the state-of-the-art predictions of semi-analytic models (SAM) of galaxy 
formation and evolution. We use the PACS Evolutionary Probe Survey (PEP) and 
Herschel Multi-tiered Extragalactic Survey (HerMES) data-sets in the COSMOS and 
GOODS-South fields, combined with broad-band photometry from UV to sub-mm, 
to obtain total (IR+UV) instantaneous star formation rates (SFRs) for individual 
Herschel galaxies up to 4, subtracted of possible active galactic nucleus (AGN) 
contamination. The comparison with model predictions shows that SAMs broadly 
reproduce the observed SFRFs up to z~2, when the observational errors on the SFR 
are taken into account. However, all the models seem to under-predict the bright- 
end of the SFRF at z> 2. The cause of this underprediction could lie in an improper 
modelling of several model ingredients, like too strong (AGN or stellar) feedback in 
the brighter objects or too low fail-back of gas, caused by weak feedback and outflows 
at earlier epochs. 

Key words: cosmology: observations - galaxies: evolution - galaxies: formation 
galaxies: luminosity function - galaxies: star-formation - infrared: galaxies. 


1 INTRODUCTION 

The study of how the star formation rate (SFR) in galaxies 
evolves with redshift provides important constraints to the 
galaxy formation and evolution theories. In particular, semi- 
analytic models (SAMs; e.g. White & Frenk 1991; Kauff- 
mann et al. 1993; Springel et al. 2001; Monaco et al. 2007; 
Guo et al. 2011; Benson 2012; Menci et al. 2012; Henriques 
et al. 2013), need to be directly compared with observa¬ 
tions to obtain insight of the relevant physical processes. The 
first and most popular SAMs are three, commonly named 
“Munich” (starting with the models of Kauffmann, White 

* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and 
with important participation from NASA 
f E-mail: carlotta.gruppioni@oabo.inaf.it 


& Guiderdoni 1993), “Durham” (beginning with the mod¬ 
els of Cole et al. 1994), and “Santa Cruz” (beginning with 
the models of Somerville & Primack 1999); more recent 
SAMs include, i.e., Croton et al. 2006; Bower et al. 2006; 
Somerville et al. 2008; Fontanot et al. 2009; Guo et al. 2011; 
Somerville et al. 2012. The main differences between these 
models lie in the prescriptions adopted for some of the most 
basic baryonic processes, such as star formation, gas cooling 
and feedback. One of the processes that must be modelled 
and compared to data is the evolution of the SFR over the 
cosmic time. However, the derivation of an accurate SFR 
from observational data is difficult, due to the many uncer¬ 
tainties involved in its reconstruction. An important source 
of uncertainty comes from dust extinction. The rest-frame 
ultraviolet (UV) fight emitted by young and massive stars, 
strictly connected to the instantaneous SFR in galaxies, is 
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strongly absorbed by dust, and re-radiated in the infrared 
(IR) bands. Dust attenuation, as well as other galaxy phys¬ 
ical properties, evolves with cosmic time and shows a peak 
between z~ 1 and 2 (e.g. Burgarella et al. 2013). Knowing 
how dust attenuation evolves in redshift is therefore crucial 
to study the redshift evolution of the SFR: to this purpose, 
combining UV information with direct observations in the 
IR region is probably the best tool to account for the total 


SFR (e.g., Kennicutt 19981. In fact, IR surveys covering a 


wide range of redshifts are extremely useful to estimate the 
global IR luminosity, since they provide a direct measure¬ 
ment of the amount of energy absorbed and re-emitted by 
dust (e.g., what is missed by UV surveys). 


Herschel (Pilbratt et al. 20101, with its 3.5-m mirror, 


has been the first telescope which allowed us to detect the 
far-IR. population to high redshifts (z~3-4) and to derive 


its rate of evolution through a detailed LF analysis (Grup- 
pioni et al.||2013 |Magnelli et al. |[2013 1 thanks to the extra- 


galactic surveys provided by the Photodetector Array Cam¬ 
era & Spectrometer ( PACS; Pogli tsch et al.|2010| and Spec¬ 
tral and Photometric Imaging Receiver (SPIRE; |Griffin et 


tionary Probe, PEP, Lutz et al. 

2011. Herschel Multi-tiered 

Extragalactic Survey, HerMES, 

Oliver et al. 2012; GOODS- 


Herschel, |Elbaz et al.||20lT| Herschel- ATLAS, | Pales et al.| 


2010). PEP and HerMES are the major Herschel Guaran¬ 


teed Time extragalactic key-projects, designed specifically 
to determine the cosmic evolution of dusty star formation 
and of the IR LF, and include the most popular and widely 
studied extragalactic fields with extensive multi-wavelength 
coverage available (deep optical, near-IR and Spitzer imag¬ 
ing and spectroscopic and photometric redshifts ; see [Berta] 
|et al.|2011| |Lutz et al.|201l||Qliver et al.|20l2| for a detailed 
description of the fields and observations). The far-IR do¬ 
main in galaxies, although potentially contaminated by the 
presence of an AGN, has been probed to be dominated by 
star formation (i.e. |Hatziminaoglou et al.|[2010| |Delvecchio| 
et al.|20l4 1 . Therefore, PEP and HerMES, and all the ancil¬ 


lary data available in the fields, give us the opportunity to 
disentangle star formation from AGN contribution and to 
study in detail the evolution of the SFR with cosmic time 
since the Universe was about a billion years old. 

In a recent paper from Pozzi et al. (2015), the observed 
IR PEP/HerMES LFs have been reproduced by means of 
a phenomenological model considering two galaxy popula¬ 
tions characterised by different evolutions, i.e. a populations 
of late-type sources and a populations of proto-spheroids. In 
the model, the IR luminosity functions (linked to the SFR) 
have been reproduced, as well as the literature K-band lumi¬ 
nosity functions (directy linked to the stellar mass), showing 
that most of the PEP-selected sources observed at z> 2 can 
be explained as progenitors of local spheroids caught during 
their formation. 

Similar evolutionary rates have been found by |Gruppioni et| 
al. (2013]l, deriving the far- and total IR (i.e. rest-frame 8- 
1000 pm) LFs from the Herschel data obtained within the 
PEP and HerMES projects up to z~ 4. Since a large fraction 
of Herschel selected objects have been found to contain an 


AGN (Berta et al.|2013 Gruppioni et al. 2013 Delvecchio et 


al. 20141, an accurate quantification of the AGN contribu¬ 


tion to the IR luminosity is needed in order to derive reliable 
SFRs (e.g., not contaminated by AGN activity) from these 


sources. Delvecchio et al. (20141, through a detailed SED 
decomposition analysis (see Berta et al.||20l 3 L have disen¬ 
tangled the contribution to the total IR luminosity due to 
AGN activity and that due to SF for the whole PEP sample. 
By starting from the work of Delvecchio et al. ( 2014[ ), but 
considering the contribution of SF only, in the present paper 
we focus on the determination of the SFR function (SFRF) 
and SFR density (SFRD) and compare the results with the 
predictions of state-of-the-art SAMs of galaxy formation and 
evolution. 

This paper is organised as follows. We discuss the PEP 
multi-wavelength catalogue in Section 2 and the theoretical 
predictions and comparison with data in Section 3; finally 
we present our conclusions in Section 4. 

Throughout this paper, we use a Chabrier (2003) initial mass 
function (IMF) and we assume a ACDM cosmology with 
H 0 = 71 km s _1 Mpc -1 , D m =0.27, and Da = 0.73 for data 
derivations. Note that, although not affecting the results of 


this paper (see, e.g., Wang et al. 20081, the considered SAMs 


use slightly different cosmologies. 


2 THE SFR FUNCTION OF HERSCHEL 
SELECTED GALAXIES 

2.1 The data-set 

We have considered the Herschel PACS (70, 100 and 160 
pm) and SPIRE (250, 350 and 500 pm) data in the COS¬ 
MOS and GOODS-S fields from the PEP and HerMES Sur¬ 
veys and all the multi-wavelength data-set associated to the 
far-IR sources. The reference sample is the PEP blind cat¬ 
alogue selected at 160 pm to the 3cr level, which consists of 
4118 and 492 sources respectively in COSMOS (to 10.2 mjy 
in ~2 deg 2 ) and GOODS-S (to 1.2 mjy in ~196 arcmin 2 . 
As described in detail by |Berta et al. (20111 and Gruppi- 


|oni et ah[ ( |2013[ > , our sources have been associated to the 
ancillary catalogues by means of a multi-band likelihood 


ratio technique (e.g., Sutherland & Saunders 19921, start¬ 
ing from the longest available wavelength (160 pm, PACS) 
and progressively matching 100pm (PACS), 70pm (PACS, 
GOODS-S only) and 24 pm (Spiteer/MIPS). In the GOODS- 
S Held, we have associated to our PEP sources the 24-pm 


catalogue by Magnelli et al. (2009) (extracted with IRAC 
3.6-pm positions as priors), that we have matched with the 
optical+near-IR+IRAC MUSIC catalogue of Grazian ct al. 
(20061, revised by Santini et al. (20091, which includes spec¬ 


troscopic and photometric redshifts. In COSMOS, we have 
matched our catalogue with the deep 24-pm sample of |Le| 


Floc’h et al. 


(20091) and with the IR AC-based catalogue of 


Ilbert et al. 820101, including optical and near-IR photome¬ 


try and photometric redshifts. In HerMES a prior source ex¬ 
traction was performed using the method presented in Rose- 
boom et al. (20111, based on MIPS-24pm positions. The 24- 
pm sources used as priors for SPIRE source extraction are 
the same as those associated with our PEP sources through 
the likelihood ratio technique. We have therefore associated 
the HerMES sources with the PEP sources by means of the 
24-pm sources matched to both samples. For most of our 
PEP sources (~87 per cent) we found a >3<r SPIRE coun¬ 
terpart in the HerMES catalogues. Redshifts (either spec¬ 
troscopic or photometric) are available for all the sources in 
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Figure 1. IR+UV SFRF estimated through the 1/V m ax method in 6 representative redshift bins, by combining the data from the 
PEP GOODS-S and COSMOS fields using the |Avni fe Bahcallj (l~980[ ) technique (black filled circles). The black solid line represents 
our best fit to our data with a modified Schechter function, while the yellow solid line is the total IR LF (without excluding AGN 
contribution through SED decomposition) obtained by Gruppioni et al. (2013), converted to a SFRF. The SFRFs of 24-/rm sources with 
log(M/M 0 )>lO in the GOODS-S by |Fontanot et al.| i 2012| | are plotted for comparison as black open squares. The SAMS predictions 
are shown as purple dotted (R-SAM), sea-green short-dashed (MORGANA), blue dot-dashed ( |De Lucia &: Bl aizot 2007J and deep-pink 
long-dashed jHenriques et a.l.|2015|) coloured lines. 


GOODS-S and for 93% of the COSMOS sample (references 
and details also in |Berta et al.||20ll] and [Gruppioni et al.| 


20131. 


2.2 The SFR Function and SFR Density 


Gruppioni et al. ( 2013|) derived the far- and total IR (i.e. 


rest-frame 8-1000 pm) LFs from the Herschel data obtained 
within the PEP and HerMES projects up to z~4. To com¬ 
pute the SFRF, we have used the same data and method 
used by [Gruppioni et al. (j2013| for deriving the total IR 
LF, but here we have subtracted - for each source individu¬ 
ally - the AGN contribution from each SED, as estimated by 
Delvecchio et al. (i20l4|, to obtain the IR luminosity due to 
SF only. In order to disentangle the possible AGN contribu¬ 
tion from that related to the host galaxy, Dclvccchio ct al. 
(|2014[) have performed a broad-band SED decomposition of 


our PEP sources using the MAGPHYS code (da Cunha et 
al.|2008 1, which is a public code using physically motivated 


templates to reproduce the observed galaxy SEDs, as mod¬ 


ified by Berta et al. (20131 to include also the AGN compo¬ 
nent (from Fritz et al.||2006 and Feltre et al.||2012 models). 


Delvecchio et al. (20141 found significant (at 99 per cent) 


contribution from AGN in 37 per cent of the PEP sources 
and used the IR luminosity of the AGN component to de¬ 
rive the AGN bolometric LF and the SMBH growth rate 
across cosmic time up to z~3. Note that |Delvecchio et alT| 


(2014) choose to derive the SMBH accretion function only 
up to 3, based on the fact that too few BH accretion data 
points were available at z>3 in order to provide an accept¬ 
able fit. On the contrary, in this work, we subtract the AGN 
component from the total IR luminosity of the sources with 
a significant AGN activity to obtain the contribution due to 
SF only (Lir). In contrast with the work of Delvecchio et 
al. (2014|, we can estimate the SFRF also in the 3<z<4.2 


interval, since the SFR data above the completeness limit 
allowed us to obtain an acceptable fit (though with larger 
uncertainties than at lower z). We have used the same cali¬ 
bration of Santini et al. (2009) and Papovich et al. (2007) to 
estimate the total instantaneous SFR (then used to derive 
the SFRF): 

FFRiR+uv/Mgyr 1 = 1.8 x 10 10 x (2.2 x Luv + Lir) (1) 


with Luv=1-5xL 27oo ^ computed from the best-fit template 
SED, and Lir =l| r. To derive the SFRF, we have used the 
1 /V m ax method ( Schmidt et al.||1968 l, combining the data 
in the two fields following Avni & Bahcall (19801. According 


to this method, the SFRF value and its uncertainty in each 
SFR bin have been computed as: 


<S>{SFR,z) = 


1 


A SFR 


:e 


Wi x 14, 


/£ 


1 


(Wi X I4iax,i) 


l](2) 
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where Knax.i is the comoving volume over which the i-tli 
galaxy could be observed, A SFR is the size of the SFR bin 
(in logarithmic scale), and Wi is the completeness correc¬ 
tion factor of the i-th galaxy. These completeness correction 
factors are a combination of the completeness corrections 
given by Berta et al. (20111, derived as described in Lutz et 


al. (20111, to be applied to each source as a function of its 


flux density, together with a correction for redsliift incom¬ 
pleteness. Additional details are given in Gruppioni et al. | 
( 2013] ). 

The resulting SFRFs in different redshift intervals are 
shown in Fig. [I] (black solid circles) and presented in Table[l] 
(together with their associated 1-a uncertainties). For com¬ 
parison, the SFRF of sources with log(M/MQ)>10 in the 
GOODS-S derived by Fontanot et al. (2012) have been also 
plotted, although the two samples have different selections 
(i.e. the Fontanot et al. 2012 sample is selected at 24 /im 
and complete in mass, while ours is flux-limited at 160-/im), 
therefore are not directly comparable at the fainter SFRs 
(affected by the sample cuts). Our best-fit solution with a 
modified Schechter function (Saunders et al. 19901 has also 
been reported (black solid line) and compared with the best- 
fit to the total IR LF of |Gruppioni et ah (|2013|) (y ellow solid 
line), converted to SFRF through the Kennicutt (19981 re¬ 
lation scaled to the Chabrier IMF (although containing also 
the AGN contribution). 

Note that, apart from the two lower redshift bins 
(2<0.45), the faint-end of our SFRFs is not constrained by 
data, therefore we can only derive the value of a at low- 2 , 
then fix it and keep that value also in the higher redshift 
bins. This Hobson’s choice implies the assumption that the 
faint-end slope does not vary with redshift. However, we note 
that at 2<0.45, where a is constrained by data, the SFRF 
obtained from the conversion of the total IR LF of Gruppioni 
et al. ( 2013| ) is flatter at low SFRs than the (IR+UV) SFRF 


(e.g. the yellow uncertainty area is lower than the faintest 
SFR data point). This can be interpreted as due either to 
a major contribution of the UV to the faint-end (almost 
negligible at higher SFRs, dominated by the IR) or/and to 
the fact that some sources might move to fainter luminosity 
bins when the AGN contribution is removed, thus steep¬ 
ening the SFRF. The latter explanation is also consistent 
with our finding that in all the redshift bins but the highest 
one (which is also the most uncertain, due to the high frac¬ 
tion of photometric redshifts), the bright-end of the SFRF 
(black solid line) is always steeper than that of the total IR 
LF converted to SFRF (yellow solid line). Since the AGN- 
dominated sources contribute mainly to the bright-end of 
the IR LF (see Gruppioni et al. 2013|, this difference is due 
to the subtraction of the AGN component from their SEDs. 

The SFRF of 24-/im selected sources with M>10 10 M@ 
by Fontanot et al. (20121 is in good agreement with our 
determination, in the common redshift and SFR range, with 
the Fontanot et al. (2012) one being flatter in the lower SFR 
common bin, likely due to the mass cut in their sample. 

By integrating the best-fitting modified Schechter func¬ 
tion to our IR+UV (IR) SFRFs down to log\o(SFR)=— 1.5 
in the different redshift bins, from 2 ~0 to 2 ~4, we have de¬ 
rived the comoving IR+UV (IR) SFR density (SFRD), as 
presented in the 10 th (last) column of Table 111 and shown 
in Fig [2] as a grey-filled (orange line-filled) area. For com¬ 
parison, other derivations from different bands are shown 
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Figure 2. Redshift evolution of the comoving SFRD. The re¬ 
sults of integrating the best-fitting curve for our observed IR+UV 
SFRD in each z-bin is shown as grey filled area (±lcr uncertainty 
locus). Top panel: Our IR+UV SFRD estimate is compared with 
other derivations: the yellow line-filled area shows the total IR Lu¬ 
minosity Density resulting from integrating the [Gruppioni et al.| 
( |2013[ > best-fitting curve (i.e. converted to SFRD without exclud¬ 
ing the AGN contribution); the orange line-filled area shows the 
uncertainty region of the only-IR SFRD from our data (after AGN 
removal); the blue line-filled area is the uncorrected for extinction 
UV SFRD by[Cucciati et al.| ( |2012[ >; the dark green filled area is 
the |Burgarella et al.| ( |2013[ ) estimate from Herschel data; the pale 
blue filled area represents the fit to optical/UV data by |Behroozi| 
|et~al. [ ( J2013| >; the magenta line is the the best-fitting function to 
IR and UV comoving SFRD by Madau & Dickinson (2014). The 
honey dew star shows the SFRD value obtained by Marchetti et 
al. (2015) from the HerMES local wide area sample, by combining 
the IR and the UV SFRs (but without excluding the AGN con¬ 
tribution). Middle panel: Our IR+UV SFRD estimate shown in 
the previous panel is compared with the SAMs predictions (same 
colours and line-styles as in the previous figures) integrated over 
the same range of SFRs (from logio(SFR)=— 1.5). Bottom panel: 
same as in the previous panel, but with our best-fitting function 
and SAMs integrated only over the SFR range covered by our 
data. 
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Table 1. PEP SFRF and SFRD 






log(«3>/Mpc 

3 dex -1 ) 




SFRDi K+xjv 
(M© yr -1 

SFRD 1K 

Mpc -3 ) 

z 




log(5FH IR _|_ uv /MQ yr - 1 ) 







-0.54-0.0 

0.04-0.5 

0.54-1.0 

1.04-1.5 

1.54-2.0 

2.04-2.5 

2.54-3.0 

3.04-3.5 



0.0—0.3 

-1.78+0.23 

-2.26+0.05 

-2.53+0.04 

-3.29+0.05 

-4.92+0.31 

-5.22+0.43 



0.025+0.005 

0.022+0.016 

0.3-0.45 


-2.52+0.14 

-2.36+0.11 

-2.85+0.08 

-4.29+0.11 

-5.22+0.31 



0.035+0.010 

0.028+0.012 

0.45-0.6 



-2.70+0.09 

-2.94+0.06 

-3.79+0.05 

-5.45+0.31 



0.049+0.014 

0.038+0.013 

0.6—0.8 



-2.38+0.19 

-2.67+0.06 

-3.51+0.03 

-4.96+0.13 



0.056+0.013 

0.039+0.015 

o 

1 

00 

6 




-3.01+0.08 

-3.26+0.04 

-4.40+0.06 

-6.17+0.43 


0.064+0.016 

0.050+0.019 

1.0-1.2 




-2.95+0.11 

-3.17+0.08 

-4.13+0.04 



0.062+0.014 

0.056+0.026 

1.2-1.7 





-2.76+0.13 

-3.90+0.04 

-5.23+0.08 


0.082+0.021 

0.051+0.025 

1.7-2.0 






-3.96+0.10 

-4.61+0.05 

-6.56+0.43 

0.071+0.019 

0.062+0.023 

2.0-2.5 






-3.54+0.13 

-4.46+0.06 

-6.09+0.19 

0.062+0.021 

0.058+0.015 

2.5-3.0 






-4.00+0.34 

-4.23+0.13 

-5.32+0.10 

0.056+0.020 

0.053+0.016 

3.0—4.2 






-4.67+0.28 

-4.94+0.20 

-5.34+0.32 

0.028+0.012 

0.016+0.012 


(i.e. the integration of the best-fitting Schechter function to 
the total IR LF, containing AGN, converted to SFR as yel¬ 
low line-filled area; the IR+UV SFRD by |Burgarella et al.| 
2013|as dark-green filled area; the fit to optical/UV data by 


Behroozi et al.||2013| as pale blue filled-area; the UV SFRD 


- uncorrected for extinction - derivation by |Cucciati et al.| 


2012 as blue line-filled area). 


From the comparison between our SFRD and previous 
derivations, we notice that the IR+UV SFRD estimated in 
this work is higher at low redshift (i.e. z<0.5), while it agrees 
within the uncertainties with the other estimates at higher 
z. The IR-only SFRD, given the larger uncertainties, is con¬ 


sistent with the optical SFRD by Behroozi et al. (20131 and 


with the Gruppioni et al. (2013) IR LD at low-z, although 


the average value is also higher than the latter estimates. 
Therefore, the low redshift difference is likely due to the 
UV SFR contribution (that in percentage is higher at low- 
z), but mostly to the AGN contribution subtraction. The 
(uncorrected for dust extinction) UV SFRD by |Cucciati et| 
|al.| ([2012]) is significantly lower than our IR or IR+UV one 
over the 0<z<3 range, while at z>3 it becomes compara¬ 
ble. This is consistent with the peak of dust extinction being 
around z~1.5-2, with dust attenuation rapidly decreasing at 
higher redshifts (>3-4; e.g., Burgarclla et al. 2013]) . In the 
top panel of Fig. [2] we also show the redshift evolution of 
the total SFRDs as obtained (from the same data sample) 
by Burgarella et al. (2013; dark-green shaded region). The 
differences (i.e. the Burgarclla et al. 2013 estimate is lower 
than the current one at z<l and slightly higher, but still 
within the uncertainties, at z>2.5) can be ascribed to the 
fact that in the previous analysis an average AGN contri¬ 
bution for each population had been subtracted from the 
total IR luminosity density (then converted to SFRD and 
summed to the UV SFRD), while in this work an accurate 
object-by-object subtraction of the IR luminosity contribu¬ 
tion due to the AGN has been performed, thanks to the 
detailed SED-fitting and decomposition of Dclvccchi o et al.| 
12014), and a proper SFRF has then been calculated from 


the obtained IR+UV SFRs. 

Finally, we note that, while the best-fitting function to 
the comoving SFRD from IR and UV data by Madau & 
Dickinson (2014) is in good agreement with our derivation 
(although slightly lower at z<0.5 and higher at z>3), at 
z>0.5 the average|Behroozi~et al. (20131 estimate is always 
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Figure 3. Total (IR+UV) SFRF (same as in Fig. black filled 
circles) in two representative redshift bins compared to SAMs 
predictions not convolved with an error of 0.3 dex (same colours 
and line-styles as in Fig. [I]). 
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higher (although consistent within the large uncertainty re¬ 
gion), maybe due to large extinction corrections. 

Finally, in Fig. [2]we show the SFRD value obtained by 
Marchetti et al. (2015) from the HerMES wide-area sample, 
by combining the SFR from IR and UV (in analogy with the 
present work), but without excluding the AGN contribution. 
The value is in good agreement with the result of|Gruppioni| 


et al. (2013) (and the others from the literature) at the same 


redshift and only marginally consistent with out IR+UV 
SFRD result. 
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3 SEMI-ANALYTICAL MODEL PREDICTIONS 

In this section, we compare the observed SFRFs with re¬ 
sults obtained with four SAMs of galaxy formation. The 
SAMs considered in this work are: the MOdel for the Rise 


of Galaxies aNd Agns (MORGANA, Monaco, Fontanot & 
Taffoni|2 007 sea-green short dashed lines in the figures), R- 
SAM (Menci et al. 2012, 2014; purple dotted lines) and two 
versions of the Munich modeQby De Lucia & Blaizot (2007; 
dot-dashed blue lines) and Henriques et al. (2015; deep-pink 
long dashed lines). Note that these SAMs use slightly differ¬ 
ent cosmologies (see the relative papers for details), although 
this does not affect the results discussed in this paper (as 
shown by Wang et al.| [2008). Moreover, all the SAMs’ pre¬ 
dictions (except the R-SAM’s ones) are mass-limited, with 
a cut at M=10 9 M0, though this selection affects only the 
low-SFRs, typically lower than those reached by our data. 
Since these models have not been "tuned” to reproduce the 
SFRFs, the results shown in this work can be considered as 
genuine SAMs predictions. 

SAMs treat the physical processes involving baryons 
(thermal state and infall/outflow of gas, star formation, feed¬ 
back, accretion onto a black hole, AGN feedback) within the 
backbone of dark matter (DM) halo merging histories, pro¬ 
duced by the gravitational collapse of DM (see Benson 2001 
for a review). We refer the reader to the papers cited above 
for all details on the models. 

Together with the observational SFRF, Fig. 1 reports, 
at the six chosen redshift ranges, the SFRFs predicted by 
the four SAMs. To take into account the errors in the ob¬ 
servational determination of SFRs, model predictions have 
been convolved with a fiducial error of 0.3 dex (assuming 
a log-normal error distribution for the SFR with an am¬ 
plitude of 0.3 dex). As discussed by Fontanot et al. (2009; 
2012), this value is roughly equal to the median formal er¬ 
ror of SFRs in the GOODS-MUSIC catalogue (Santini et 
al. 2009) and allows us to determine the gross effect of (ran¬ 
dom) uncertainties in SFR determinations. We consider this 
value suitable also for our far-IR-based SFR. 

At low redshift (*<1.7, the three upper panels), all but 
MORGANA model predictions, are remarkably similar on 
the bright end. The MORGANA model at 2 <0.8 shows an 
excess of very bright sources, that is connected to the spe¬ 
cific model of radio-mode AGN feedback, where accretion 
onto the central black hole takes place from the cold, star 
forming gas in the bulge, so the suppression of star forma¬ 
tion is partial. All models tend to give a power-law tail for 
the bright-end of the SFRF, that is broadly compatible with 
the analytical extrapolation of the observed SFRF presented 
in Section 2. The broadly good agreement between models 
and data at the bright-end breaks in the very important 
redshift range from 2 to 3, with the very steep SFRFs pro¬ 
duced by all models becoming apparently different from the 
analytic fit to data in the highest redshift bin. However, 
the slope of this analytic fit is determined by * = 0 — 0.3 
data, while at z > 2, there is clear observational evidence 
of a steep (or even steepening of the) UV luminosity func¬ 
tion (e.g. Cucciati et al. 2012). The consistency at least 
of one of these models (e.g. MORGANA) with the LF of 

1 obtained from the publicly available database 
http://gavo.mpa-garching.mpg.de/Millennium/ 


Lyman-break galaxies has been investigated by Lo Faro et 
al. (2009). At 2>2, the models underpredict the bright-end 
of the SFRF by a factor of ~0.2—0.3 dex. This result de¬ 
pends sensitively on the adopted modelling of observational 
error. Figure 3 shows the comparison of models and data for 
two redshift bins at 2^0 and 2 ~2, without convolving with 
0.3 dex error. While the comparison at 2=0 remains accept¬ 
able, the disagreement at 2 ~ 2 worsens dramatically. Only 
extreme assumptions on the error on SFRs would allow to 
recover the 2 ~2 SFRF. 

We stress that in the highest redshift bin (3<2<4.2), 
where the knee of the observed SFRF is not clearly de¬ 
tectable and the shape seems to change significantly, the 
fraction of photometric redshifts and of the power-law SED 
AGN (with an uncertain photo -2 determination) is higher 
than at lower 2 ’s. For this reason, this 2 -bin is more affected 
by uncertainties than the other ones and the discrepancy 
with model predictions here must be taken just as indica¬ 
tions (to be furtherly verified). 

The inability of SAMs to reproduce the bright-end of 
the Herschel LF at 2 > 2 had previously been found bylNiemi 


et al. (2012 I, though based on very preliminary Herschel LF 


results. It might be connected to the tendency of models to 
underestimate the main sequence of star-forming galaxies at 
2 ~2 (e.g., Dutton et al. 2010). Because SFRs are determined 
by the complex pattern of gas inflows and outflows, in and 
out DM halos, and because inflows are determined by out¬ 
flows taking place at previous times, identifying the cause 
of this disagreement is not easy. One possibility could be 
excessive feedback, possibly from AGN either in the radio 
mode (cooling and thus SFR is over-suppressed in these ob¬ 
jects) or in the quasar mode (quasar-triggered outflows limit 
SFR); however, models with very different implementations 
of AGN feedback are showing the same problem. Alterna¬ 
tively, an excessive formation of stars in low-mass galaxies 


(Fontanot et al. 2009; Lo Faro et al. 2009; Weinmann, Neis- 
|tein fe Dekel|20iT] ) could lock too much gas in stars instead 
of ejecting it from DM halos and making it available for later 
star formation. 

What is interesting to note is that, despite the very dif¬ 
ferent prescriptions for SFR and stellar feedback, all models 
provide similar predictions (at least up to 2~3) in terms 
of number density of objects with moderate SFRs. Proba¬ 
bly this reflects the fact that all the models are calibrated 
to reproduce the knee of the 2=0 galaxy luminosity/mass 
function. 

At low SFRs (< 1 Mq yr -1 ), typically below the com¬ 
pleteness limit of Herschel data, models start to separate, 
with R-SAM and Henriques et al. (2015) giving respectively 
the highest and lowest SFRFs. These differences may be due 
to a number of features: different modelling of stellar and 
AGN feedback, calibration procedure, treatment of merger 
trees (analytical or based on simulations). The turnover at 
low SFRs observed in all models but R-SAM is mostly due to 
incompleteness: models are complete in DM halo mass, that 
is tightly correlated with stellar mass, while the correlation 
with SFR is much broader and time-dependent, and this 
creates a very broad cutoff. While at 2 = 0 model SFRFs 
tend to be lower than the observed ones (as noticed, e.g., 
by Fontanot et al. 2009), at higher redshifts (but still be¬ 
low 2~1.7) they overshoot by a large factor the SFRF of 
Fontanot et al. (2012), derived for 24-/im selected sources 
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with M* > 10 10 M@ in the GOODS-S (but are consistent, 
up to 2 ~2, with the analytic fit to our SFRF). Part of the 
difference with Fontanot et al. (20121 is due to the selection 


in stellar mass done in that paper, which produce a flat¬ 
tening (double-peaked) of the faint-end, although the same 
paper shows that MORGANA cannot predict the drop in 
the SFRF at ~ 1 M@ yr -1 . 

The overprediction of models with respect to the an¬ 
alytic fit to data observed at logio(SFR)<l — 2, becom¬ 
ing more and more important with increasing redshift (al¬ 
though the faint-end slope of the observed SFRSs is fixed at 
the value found at 2 ~ 0 ), if confirmed by deeper data, could 
be due to the well known excess of low/intermediate mass 
galaxies predicted by most SAMs with respect to observed 
mass functions (see, i.e., Somerville & Primack 1999; Cole 
et al. 2000; Menci et al. 2002; Croton et al. 2006). This as¬ 
pect is a result of the small-scale power excess typical of the 
A — CDM power spectrum (e.g., Moore et al.|1999||Klypin 


et al. |1999 Menci et al.|2012| Calura et al.|2014 |. Some au¬ 

thors have suggested that a strong feedback from exploding 
supernovae can help limiting this excess, and alleviate the 
discrepancy between data and models at the smallest scales 


(e.g., Croton et al. |2006| Pontzen fe Governato|2dl2 |. 

The integrals of the SFRFs (SFRD, shown in Figure 
2, middle and bottom panels) confirm and better quantify 
the trends shown in Figure 1. The SFRDs obtained by in¬ 
tegrating the SAMs over the same luminosity range as our 
best-fit SFRFs, are marginally low at z = 0 for three models 
out of four (the MORGANA prediction is very close to the 
observed one, but this is driven by the bright excess visible 
in Fig. 1). MORGANA and the two Munich models follow 
the evolution of the observed SFRD up to 2 ~2, while the R- 
SAM prediction is significantly steeper, overpredicting the 
data estimate from 2 ~0.6 up to the higher redshifts. From 
2=2 to 2=3 all models but |Henriq ues elTal. (2015) overpre¬ 
dict the observed SFRD (at 1 -<t uncertainty level) by a factor 
>2, with R-SAM diverging by a factor as high as 7-8. Note 
that if we consider a more conservative uncertainty level for 
our data (i.e., 3 -ct) the disagreement would be much less or 
even negligible. 

Since the observed SFRD shown in the upper and mid¬ 
dle panels of Fig. [2] is obtained by integrating the SFRFs 
down to SFR values not covered by data, we have also com¬ 
puted the SFRD (from both data and models) by perform¬ 
ing the integration just in the SFR range where PEP data 
are available (as shown in the bottom panel). Here the trend 
of models underpredicting the observed SFRD at high 2 is 
even more evident, with data values starting to be underes¬ 
timated by SAMs at 2~1.5-1.8, but getting closer to model 


predictions (especially with the Henriques et al. 2015 and 
R-SAM ones) at 2~3.5 (although the better agreement of 
the 3.2<2<4 SFRD might be due to the compensation of 
the under-estimate of the bright-end and the over-estimate 
of the faint-end of the SFRF). 

We can conclude that our results at high -2 confirm a 
tension between models and data that possibly points to 
a problem related to the common assumptions that are at 
the base of galaxy formation within the ACDM cosmogony. 
Given the good agreement between models and data at low- 
2 , we can interpret these z>2 tensions as a consequence of 
model assumptions (like, e.g., the SF law and efficiency, IMF 
shape) and parameters being calibrated with local observa¬ 


tions, and then assumed invariant at higher redshift. It is in¬ 
deed possible for some of these analytical approximations to 
have an intrinsic or acquired redshift dependence (through 
the evolution of the physical properties of model galaxies, 
see, e.g., Fontanot 2014, where the predicted SFRFs for 
SAMs with variable IMFs are discussed). Nonetheless, the 
fact that the largest discrepancies between data and models 
are seen for the 2<2<3 interval, which represents a peak 
for both the cosmic SFR and the BH accretion, points to 
the treatment of SFR in extreme environments as a likely 
source of the tension. Indeed, the modelling of extreme envi¬ 
ronments, such those associated with the strongest starburst 
is still highly uncertain, with many theoretical studies sug¬ 
gesting a different star formation regime for these systems 
(see, e.g., Somerville, Primack fe Faber|200l Hopkins et al. 


2010 ). 


On the other hand, also possible source of uncertainty in the 
data (affecting mainly the highest redshift interval) might 
contribute to the tension, increasing the bright-end of the 
observed SFRF: they can be due to wrong photometric red- 
shifts, bright IR sources at confusion level with flux en¬ 
hanced by “blending”, mis-identification of lensed galaxies 
(the latter only to a very minor extent, since the 160-/rm se¬ 
lection should be much less affected by lensed galaxies than 
sub-mm wavelength ones; e.g. Negrello et al. [2007). 

In the future, deeper far-IR observations will be funda¬ 
mental to explore the faintest end of the SFRF, whereas a 
study of the bright end at larger redshifts will provide tighter 
constraints on the feedback processes regulating star forma¬ 
tion in the brightest galaxies. A combined study of SFRF 
and mass function will be also crucial in order to fully un¬ 
derstand and tune the single processes considered in SAMs. 


4 CONCLUSIONS 


Starting from the far-IR PEP data considered by Gruppi- 
oni et al. (2013), in this paper we investigate the evolution 
of the SFRF in the redshift range 0.1<2<4 and compare it 
with theoretical results from various semi-analytical mod¬ 


els of galaxy formation (MORGANA, R-SAM, De Lucia & 
Blaizot 2007 and Henriques et al. I 20151. To compute the 


SFRF, we have subtracted the AGN contribution estimated 


by Delvecchio et al. (20141 from each SED, to obtain the 


IR luminosity due to SF only. Then, we have obtained the 
total instantaneous SFR by combining the IR SF luminosity 
with the UV derived one, to obtain the total (IR+UV) SFR, 
which we have considered to derive the SFRF. 

The conclusions of our work can be summarised as follows. 


• We find a generally good agreement between the ob¬ 
served and predicted SFRFs up to 2 ~2 (once the obser¬ 
vational errors are taken into account in SAMs), with the 
exception of MORGANA, showing a high-SFR excess at 
2 <0.8. This result implies that theoretical models, despite 
the different prescriptions, are able to reproduce the space- 
density evolution of the IR luminous galaxies from the SFRD 
peak epoch up to now. 

• At 2 J 2, all the models start to under predict the bright- 
end of the SFRF. This can be due to improper modelling of 
several ingredients that determine the inflow/outflow pat¬ 
terns of gas in/from DM halos, like AGN feedback, limiting 
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SFR in the largest galaxies, or even inefficient feedback from 
a previous generation of galaxies. 

• Our data are able to constrain the low-SFR end of the 
(IR+UV) SFRF only at low -2 (2<0.45). In this range of 
redshift, the observed slope is consistent with model predic¬ 
tions, but is steeper than the slope of the total IR LF (con¬ 
taining also AGN contribution) by Gruppioni et al. (20131. 
However, at intermediate/high-z, Herschel data do not sam¬ 
ple the low-SFR end of the SFRF, where SAM predictions 
differ most: our data are not deep enough to allow us to 
distinguish between the different approaches to SF and stel¬ 
lar feedback considered by the different models. Additional 
sources of uncertainty affecting the models could be due to 
the fact that no evolution with redshift is considered for local 
relations and functions, as SFR laws and IMF. Finally, also 
data at high redshift could be affected by wrong photomet¬ 
ric redshifts and source confusion, contributing to enhance 
the bright-end of the SFRF, therefore the discrepancy with 
model predictions. 


In this work we have shown that SFRF may help putting 
stringent constraints on the physical processes modelled in 
SAMs, especially if extended to low-SFRs, while a study of 
the bright-end at larger redshifts will provide tighter con¬ 
straints on the feedback processes regulating star formation 
in the brightest galaxies. A combined study of SFRF and 
mass function will be crucial in order to fully understand 
and tune the single processes modelled in SAMs and to have 
a global picture of the evolution of SFR and mass growth in 
galaxies. 
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